The goal here is to use conventional alpha diversity metrics to see how Chao1 richness, shannon diversity and evenness change across samples and to compare those to the values seen using breakaway in the AlphaDiversity.Rmd file
Run AlphaDiversity in scratchnotebooks That file calculates richness in breakawy which I will combine here
#source(here::here("RScripts", "InitialProcessing_3.R"))
source(here::here("RLibraries", "ChesapeakePersonalLibrary.R"))
Registered S3 methods overwritten by 'dbplyr':
method from
print.tbl_lazy
print.tbl_sql
── Attaching packages ───────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.2 ──✔ ggplot2 3.3.6 ✔ purrr 0.3.4
✔ tibble 3.1.8 ✔ dplyr 1.0.10
✔ tidyr 1.2.1 ✔ stringr 1.4.1
✔ readr 2.1.3 ✔ forcats 0.5.2 ── Conflicts ──────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ksource(here::here("ActiveNotebooks", "BreakawayAlphaDiversity.Rmd"))
processing file: /home/jacob/Projects/ChesapeaakeMainstemAnalysis/ActiveNotebooks/BreakawayAlphaDiversity.Rmd
|
| | 0%
|
|...... | 5%
|
|........... | 10%
|
|................. | 14%
|
|....................... | 19%
|
|............................. | 24%
|
|.................................. | 29%
|
|........................................ | 33%
|
|.............................................. | 38%
|
|................................................... | 43%
|
|......................................................... | 48%
|
|............................................................... | 52%
|
|..................................................................... | 57%
|
|.......................................................................... | 62%
|
|................................................................................ | 67%
|
|...................................................................................... | 71%
|
|........................................................................................... | 76%
|
|................................................................................................. | 81%
|
|....................................................................................................... | 86%
|
|............................................................................................................. | 90%
|
|.................................................................................................................. | 95%
|
|........................................................................................................................| 100%
output file: /tmp/RtmpIdhQyV/file353597baadc
Registered S3 method overwritten by 'data.table':
method from
print.data.table
Registered S3 methods overwritten by 'htmltools':
method from
print.html tools:rstudio
print.shiny.tag tools:rstudio
print.shiny.tag.list tools:rstudio
Attaching package: ‘flextable’
The following object is masked from ‘package:purrr’:
compose
Attaching package: ‘ftExtra’
The following object is masked from ‘package:flextable’:
separate_header
Warning: Assuming taxa are rows
library(vegan)
Loading required package: permute
Loading required package: lattice
This is vegan 2.6-3
library(cowplot)
library(flextable)
library(ftExtra)
This file is dedicated to conventional, non div-net/breakaway stats, since breakaway seems to choke on this data.
Reshape back into an ASV matrix, but this time correcting for total abundance
nonSpikes %>% head
raDf <- nonSpikes %>% pivot_wider(id_cols = ID, names_from = ASV, values_from = RA, values_fill = 0)
raMat <- raDf %>% column_to_rownames("ID")
raMat1 <- raMat %>% as.matrix()
countMat <- nonSpikes %>%
pivot_wider(id_cols = ID, names_from = ASV, values_from = reads, values_fill = 0) %>%
column_to_rownames("ID") %>% as.matrix()
seqDep <- countMat %>% apply(1, sum)
min(seqDep)
[1] 852
sampleRichness <- rarefy(countMat, min(seqDep))
rarefy everything to the minimum depth (852)
countRare <- rrarefy(countMat, min(seqDep))
Gamma diversity
specpool(countRare)
Doesn’t finish
#specpool(countMat)
All richness estimates
richnessRare <- estimateR(countRare)
Shannon diversity
shan <- diversity(countRare)
shan
3-1-B-0-2 3-1-B-1-2 3-1-B-180 3-1-B-20 3-1-B-5 3-1-B-500 3-1-B-53 3-1-S-0-2 3-1-S-1-2 3-1-S-180 3-1-S-20
4.358550 5.057168 4.645505 5.866334 5.195276 3.865637 5.506752 4.519812 4.961172 4.823132 5.314521
3-1-S-5 3-2-B-0-2 3-2-B-1-2 3-2-B-180 3-2-B-20 3-2-B-5 3-2-B-500 3-2-B-53 3-2-S-0-2 3-2-S-1-2 3-2-S-180
4.872777 4.372510 4.698438 4.575696 5.040413 5.437996 5.035825 4.843247 3.830293 4.950842 4.746431
3-2-S-20 3-2-S-5 3-2-S-500 3-2-S-53 3-3-B-0-2 3-3-B-1-2 3-3-B-180 3-3-B-20 3-3-B-5 3-3-B-500 3-3-B-53
4.889092 5.101738 4.806509 4.395303 4.241859 4.955895 3.457605 5.716900 5.220893 5.415232 5.436138
3-3-S-180 3-3-S-20 3-3-S-500 3-3-S-53 4-3-B-0-2 4-3-B-1-2 4-3-B-180 4-3-B-20 4-3-B-5 4-3-B-500 4-3-B-53
5.025812 4.927400 4.933217 4.319009 4.323957 4.946950 4.417338 4.456536 4.658738 4.345363 4.207244
4-3-O-1-2 4-3-O-180 4-3-O-5 4-3-O-500 4-3-O-53 4-3-S-0-2 4-3-S-180 4-3-S-20 4-3-S-500 4-3-S-53 5-1-S-1-2
5.033048 4.605692 5.147523 3.933170 4.501961 2.805419 4.558078 4.656982 4.669918 4.517247 4.362809
5-1-S-180 5-1-S-20 5-1-S-5 5-1-S-500 5-1-S-53 5-5-B-0-2 5-5-B-180 5-5-B-5 5-5-B-500 5-5-B-53 5-5-S-180
4.545527 4.163814 3.991503 4.054552 4.042100 4.632035 5.108552 5.484409 5.154378 4.903917 4.280018
5-5-S-5 5-5-S-500 5-5-S-53 C_5P1B_0P2 C_5P1B_180 C_5P1B_1P2 C_5P1B_20 C_5P1B_500 C_5P1B_53
4.934133 4.921967 4.306067 4.286760 4.825743 4.766759 5.454785 4.732364 5.196367
Evenness
pielouJ <- shan/richnessRare["S.chao1",]
pielouJ
3-1-B-0-2 3-1-B-1-2 3-1-B-180 3-1-B-20 3-1-B-5 3-1-B-500 3-1-B-53 3-1-S-0-2 3-1-S-1-2 3-1-S-180
0.012650735 0.008180042 0.010628825 0.003710882 0.006572407 0.042372070 0.008906998 0.010251663 0.006769361 0.008460575
3-1-S-20 3-1-S-5 3-2-B-0-2 3-2-B-1-2 3-2-B-180 3-2-B-20 3-2-B-5 3-2-B-500 3-2-B-53 3-2-S-0-2
0.007758092 0.008185653 0.013831141 0.010000060 0.012731121 0.006334166 0.006947808 0.007370699 0.010310991 0.017040735
3-2-S-1-2 3-2-S-180 3-2-S-20 3-2-S-5 3-2-S-500 3-2-S-53 3-3-B-0-2 3-3-B-1-2 3-3-B-180 3-3-B-20
0.008120976 0.010779503 0.007686427 0.009669012 0.008608205 0.010571985 0.016670593 0.009355428 0.062019818 0.004974462
3-3-B-5 3-3-B-500 3-3-B-53 3-3-S-180 3-3-S-20 3-3-S-500 3-3-S-53 4-3-B-0-2 4-3-B-1-2 4-3-B-180
0.007297796 0.006684979 0.005663804 0.007174607 0.009110878 0.011223098 0.013101420 0.009129840 0.007506505 0.009866906
4-3-B-20 4-3-B-5 4-3-B-500 4-3-B-53 4-3-O-1-2 4-3-O-180 4-3-O-5 4-3-O-500 4-3-O-53 4-3-S-0-2
0.009604603 0.009595357 0.009468610 0.009368606 0.008100300 0.012180624 0.008613659 0.015215359 0.009862989 0.124685269
4-3-S-180 4-3-S-20 4-3-S-500 4-3-S-53 5-1-S-1-2 5-1-S-180 5-1-S-20 5-1-S-5 5-1-S-500 5-1-S-53
0.010737086 0.007742281 0.012267734 0.011180515 0.010146067 0.011378040 0.015622114 0.016866247 0.011078011 0.014276923
5-5-B-0-2 5-5-B-180 5-5-B-5 5-5-B-500 5-5-B-53 5-5-S-180 5-5-S-5 5-5-S-500 5-5-S-53 C_5P1B_0P2
0.008126376 0.008491069 0.006870717 0.009492710 0.007819712 0.013804485 0.011750911 0.009099696 0.016369217 0.006173593
C_5P1B_180 C_5P1B_1P2 C_5P1B_20 C_5P1B_500 C_5P1B_53
0.007335071 0.007715700 0.003990014 0.007887274 0.004981453
diversityData <- sampleData %>%
left_join(richnessRare %>% t() %>% as.data.frame() %>% rownames_to_column("ID"), by = "ID") %>%
left_join(shan %>% enframe(name = "ID", value = "shannonH"), by = "ID") %>%
left_join(pielouJ %>% enframe(name = "ID", value = "pielouJ"), by = "ID") %>%
arrange(Size_Class)
Parameters for all plots
plotSpecs <- list(
facet_wrap(~Depth, ncol = 1) ,
theme_bw(base_size = 16) ,
geom_point(size = 4) ,
geom_path(aes(color = as.factor(Station))) ,
scale_x_log10(breaks = my_sizes, labels = as.character(my_sizes)) ,
#scale_y_log10nice() ,
scale_shape_manual(values = rep(21:25, 2)) ,
scale_fill_viridis_d(option = "plasma") ,
scale_color_viridis_d(option = "plasma") ,
labs(x = expression(paste("Particle Size (", mu, "m)"))) ,
theme(legend.position = "none",
plot.margin = unit(c(0, 0, 0, 0), "cm"),
axis.title.x = element_blank(),
axis.text.x = element_text(angle = 90, vjust = .5),
axis.title.y = element_text(margin = unit(c(3, 3, 3, 3), "mm"), vjust = 0))
)
Observed species counts, on rarefied data
plotObs <- diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = S.obs, shape = as.factor(Station), fill = as.factor(Station))) +
plotSpecs +ylab("Observed ASVs (Rarefied)")#+ scale_y_log10()
plotObs
Estemated Chao1 Richness
plotChao1 <- diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = S.chao1, shape = as.factor(Station), fill = as.factor(Station))) +
plotSpecs +
geom_errorbar(aes(ymin = S.chao1 -2 * se.chao1, ymax = S.chao1 + 2* se.chao1), width = -.1) +
scale_y_log10() +
ylab("Richness (Chao1)")
plotChao1
Shannon diversity
plotShan <- diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = shannonH, shape = as.factor(Station), fill = as.factor(Station))) +
plotSpecs +
ylab("Diversity (Shannon H)") +
lims(y = c(2.5, 6))
plotShan
Evenness
plotPielou <- diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = pielouJ, shape = as.factor(Station), fill = as.factor(Station))) +
plotSpecs +scale_y_log10() +ylab("Evenness (PielouJ)")
plotPielou
All plots together
plotAlpha <- plot_grid(plotObs, plotChao1, plotShan, plotPielou, nrow = 1, labels = LETTERS)
plotAlpha
ggsave(here::here("Figures", "ConventionalAlpha.png"), plotAlpha, width = 11, height = 4)
Rarefied observed species numbers
obsMod <- lm(S.obs ~ log(Size_Class) + I(log(Size_Class)^2) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
summary(obsMod)
Call:
lm(formula = S.obs ~ log(Size_Class) + I(log(Size_Class)^2) +
I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
Residuals:
Min 1Q Median 3Q Max
-223.262 -37.787 -0.657 47.789 195.661
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 267625.988 133744.629 2.001 0.049325 *
log(Size_Class) 32.892 8.103 4.059 0.000128 ***
I(log(Size_Class)^2) -6.311 1.508 -4.186 8.24e-05 ***
lat -13937.556 6969.260 -2.000 0.049453 *
I(lat^2) 181.530 90.726 2.001 0.049344 *
depth 4.420 3.209 1.377 0.172926
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 76.18 on 69 degrees of freedom
Multiple R-squared: 0.2636, Adjusted R-squared: 0.2103
F-statistic: 4.94 on 5 and 69 DF, p-value: 0.0006353
Rarified chao1 estimates
chao1Mod <- lm(S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
summary(chao1Mod)
Call:
lm(formula = S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2) +
I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
Residuals:
Min 1Q Median 3Q Max
-527.97 -117.95 -33.57 128.61 860.90
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 842634.852 420822.680 2.002 0.04918 *
log(Size_Class) 78.301 25.496 3.071 0.00305 **
I(log(Size_Class)^2) -16.026 4.744 -3.378 0.00120 **
lat -43914.772 21928.525 -2.003 0.04915 *
I(lat^2) 572.134 285.467 2.004 0.04898 *
depth 18.160 10.098 1.798 0.07650 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 239.7 on 69 degrees of freedom
Multiple R-squared: 0.1924, Adjusted R-squared: 0.1339
F-statistic: 3.289 on 5 and 69 DF, p-value: 0.01011
As above but without latitude and depth
chao1ModSimple <- lm(S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2), data = diversityData)
summary(chao1ModSimple)
Call:
lm(formula = S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2),
data = diversityData)
Residuals:
Min 1Q Median 3Q Max
-467.36 -149.31 -38.88 129.05 936.52
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 554.763 43.318 12.807 < 2e-16 ***
log(Size_Class) 78.966 25.726 3.070 0.00302 **
I(log(Size_Class)^2) -16.380 4.785 -3.423 0.00103 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 242.2 on 72 degrees of freedom
Multiple R-squared: 0.14, Adjusted R-squared: 0.1161
F-statistic: 5.859 on 2 and 72 DF, p-value: 0.004389
shanMod <- lm(shannonH ~ log(Size_Class) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
summary(shanMod)
Call:
lm(formula = shannonH ~ log(Size_Class) + I(log(Size_Class)^2) +
lat + I(lat^2) + depth, data = diversityData)
Residuals:
Min 1Q Median 3Q Max
-1.32065 -0.18411 0.01949 0.27550 0.85829
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.667e+03 7.858e+02 2.122 0.0375 *
log(Size_Class) 2.024e-01 4.761e-02 4.250 6.57e-05 ***
I(log(Size_Class)^2) -3.673e-02 8.858e-03 -4.146 9.46e-05 ***
lat -8.661e+01 4.095e+01 -2.115 0.0380 *
I(lat^2) 1.127e+00 5.331e-01 2.115 0.0381 *
depth 1.869e-02 1.886e-02 0.991 0.3250
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.4476 on 69 degrees of freedom
Multiple R-squared: 0.3035, Adjusted R-squared: 0.2531
F-statistic: 6.014 on 5 and 69 DF, p-value: 0.0001125
pielouMod <- lm(pielouJ ~ log(Size_Class) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
summary(pielouMod)
Call:
lm(formula = pielouJ ~ log(Size_Class) + I(log(Size_Class)^2) +
lat + I(lat^2) + depth, data = diversityData)
Residuals:
Min 1Q Median 3Q Max
-0.016294 -0.004747 -0.002583 0.000916 0.099924
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.207e+01 2.632e+01 -0.839 0.4046
log(Size_Class) -3.952e-03 1.595e-03 -2.478 0.0157 *
I(log(Size_Class)^2) 6.757e-04 2.967e-04 2.277 0.0259 *
lat 1.149e+00 1.372e+00 0.838 0.4049
I(lat^2) -1.494e-02 1.786e-02 -0.837 0.4055
depth -3.002e-04 6.316e-04 -0.475 0.6361
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.01499 on 69 degrees of freedom
Multiple R-squared: 0.09657, Adjusted R-squared: 0.0311
F-statistic: 1.475 on 5 and 69 DF, p-value: 0.2092
uomisto H (2010a). “A diversity of beta diver- sities: straightening up a concept gone awry. 1. Defining beta diversity as a function of alpha and gamma diversity.” Ecography, 33, 2–2
predict(obsMod, se.fit = TRUE)
$fit
1 2 3 4 5 6 7 8 9 10 11 12 13 14
225.1538 225.1538 194.8219 194.8219 205.3509 156.7467 156.7467 184.4646 198.4497 300.2260 300.2260 269.8941 269.8941 280.4231
15 16 17 18 19 20 21 22 23 24 25 26 27 28
231.8189 231.8189 259.5367 259.5367 331.0301 331.0301 300.6982 300.6982 311.2272 262.6230 262.6230 290.3408 304.3259 304.3259
29 30 31 32 33 34 35 36 37 38 39 40 41 42
336.3394 336.3394 306.0076 306.0076 316.5366 316.5366 267.9323 267.9323 295.6502 295.6502 325.5520 295.2201 295.2201 305.7491
43 44 45 46 47 48 49 50 51 52 53 54 55 56
305.7491 257.1449 257.1449 257.1449 284.8628 284.8628 298.8479 298.8479 295.0645 295.0645 264.7326 264.7326 275.2616 275.2616
57 58 59 60 61 62 63 64 65 66 67 68 69 70
226.6574 226.6574 226.6574 254.3753 254.3753 268.3604 268.3604 255.1192 224.7873 224.7873 235.3163 235.3163 186.7121 186.7121
71 72 73 74 75
186.7121 214.4299 214.4299 228.4151 228.4151
$se.fit
[1] 26.65044 26.65044 25.92384 25.92384 25.78549 27.43275 27.43275 28.97311 33.61118 18.69492 18.69492 18.08628 18.08628 16.98512
[15] 19.78897 19.78897 21.23350 21.23350 19.01263 19.01263 18.55660 18.55660 17.02002 19.83509 19.83509 20.93929 27.82376 27.82376
[29] 19.01048 19.01048 18.52392 18.52392 16.74918 16.74918 19.36508 19.36508 20.31747 20.31747 18.32932 17.69829 17.69829 15.79767
[43] 15.79767 18.21194 18.21194 18.21194 19.19052 19.19052 26.22172 26.22172 19.07012 19.07012 18.19426 18.19426 16.47139 16.47139
[57] 18.19575 18.19575 18.19575 19.25451 19.25451 25.78762 25.78762 24.15321 23.21205 23.21205 22.06503 22.06503 22.85093 22.85093
[71] 22.85093 23.83784 23.83784 28.86054 28.86054
$df
[1] 69
$residual.scale
[1] 76.17928
diversityData$pred_obs = predict(obsMod, se.fit = TRUE)$fit
diversityData$se_obs = predict(obsMod, se.fit = TRUE)$se.fit
plotSpecs2 <- list(
facet_wrap(~Depth, ncol = 1) ,
theme_bw(base_size = 16) ,
#geom_point(size = 4) ,
geom_path(aes(color = as.factor(Station))) ,
scale_x_log10(breaks = my_sizes, labels = as.character(my_sizes)) ,
#scale_y_log10nice() ,
scale_shape_manual(values = rep(21:25, 2)) ,
scale_fill_viridis_d(option = "plasma") ,
scale_color_viridis_d(option = "plasma") ,
labs(x = expression(paste("Particle Size (", mu, "m)"))) ,
theme(legend.position = "none",
plot.margin = unit(c(0, 0, 0, 0), "cm"),
axis.title.x = element_blank(),
axis.text.x = element_text(angle = 90, vjust = .5),
axis.title.y = element_text(margin = unit(c(3, 3, 3, 3), "mm"), vjust = 0))
)
plotObs_pred <- diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = pred_obs, shape = as.factor(Station), fill = as.factor(Station))) +
geom_segment(aes(y = pred_obs - 2 * se_obs, yend = pred_obs + 2 * se_obs, xend = Size_Class, color = as.factor(Station)), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"), alpha = 0.5) +
plotSpecs2 + ylab("Predicted ASVs")
plotObs_pred
predict(chao1Mod, se.fit = TRUE)
$fit
1 2 3 4 5 6 7 8 9 10 11 12 13 14
461.6720 461.6720 360.3418 360.3418 441.7518 295.7073 295.7073 410.5479 376.5224 642.9478 642.9478 541.6177 541.6177 623.0276
15 16 17 18 19 20 21 22 23 24 25 26 27 28
476.9831 476.9831 591.8238 591.8238 713.7117 713.7117 612.3816 612.3816 693.7915 547.7470 547.7470 662.5877 628.5621 628.5621
29 30 31 32 33 34 35 36 37 38 39 40 41 42
719.9454 719.9454 618.6152 618.6152 700.0252 700.0252 553.9807 553.9807 668.8213 668.8213 687.4544 586.1243 586.1243 667.5342
43 44 45 46 47 48 49 50 51 52 53 54 55 56
667.5342 521.4897 521.4897 521.4897 636.3303 636.3303 602.3048 602.3048 603.6376 603.6376 502.3074 502.3074 583.7174 583.7174
57 58 59 60 61 62 63 64 65 66 67 68 69 70
437.6728 437.6728 437.6728 552.5135 552.5135 518.4879 518.4879 496.8539 395.5237 395.5237 476.9337 476.9337 330.8892 330.8892
71 72 73 74 75
330.8892 445.7298 445.7298 411.7043 411.7043
$se.fit
[1] 83.85465 83.85465 81.56842 81.56842 81.13312 86.31617 86.31617 91.16284 105.75637 58.82290 58.82290 56.90783
[13] 56.90783 53.44308 62.26529 62.26529 66.81044 66.81044 59.82256 59.82256 58.38766 58.38766 53.55289 62.41039
[25] 62.41039 65.88474 87.54647 87.54647 59.81580 59.81580 58.28484 58.28484 52.70071 52.70071 60.93154 60.93154
[37] 63.92818 63.92818 57.67254 55.68702 55.68702 49.70679 49.70679 57.30320 57.30320 57.30320 60.38230 60.38230
[49] 82.50570 82.50570 60.00344 60.00344 57.24758 57.24758 51.82664 51.82664 57.25226 57.25226 57.25226 60.58363
[61] 60.58363 81.13983 81.13983 75.99721 73.03587 73.03587 69.42682 69.42682 71.89962 71.89962 71.89962 75.00489
[73] 75.00489 90.80864 90.80864
$df
[1] 69
$residual.scale
[1] 239.6954
diversityData$pred_chao1 = predict(chao1Mod, se.fit = TRUE)$fit
diversityData$se_chao1 = predict(chao1Mod, se.fit = TRUE)$se.fit
plotChao1_pred <- diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = pred_chao1, shape = as.factor(Station), fill = as.factor(Station))) +
geom_segment(aes(y = pred_chao1 - 2 * se_chao1, yend = pred_chao1 + 2 * se_chao1, xend = Size_Class, color = as.factor(Station)), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"), alpha = 0.5) +
plotSpecs2 + ylab("Predictd Richness (Chao1)") + scale_y_log10()
plotChao1_pred
predict(shanMod, se.fit = TRUE)
$fit
1 2 3 4 5 6 7 8 9 10 11 12 13 14
4.500319 4.500319 4.326415 4.326415 4.297065 3.993768 3.993768 4.116157 4.383554 4.956811 4.956811 4.782907 4.782907 4.753557
15 16 17 18 19 20 21 22 23 24 25 26 27 28
4.450260 4.450260 4.572649 4.572649 5.151681 5.151681 4.977777 4.977777 4.948427 4.645130 4.645130 4.767519 5.034915 5.034915
29 30 31 32 33 34 35 36 37 38 39 40 41 42
5.197726 5.197726 5.023822 5.023822 4.994472 4.994472 4.691175 4.691175 4.813564 4.813564 5.145591 4.971687 4.971687 4.942337
43 44 45 46 47 48 49 50 51 52 53 54 55 56
4.942337 4.639040 4.639040 4.639040 4.761429 4.761429 5.028825 5.028825 4.981512 4.981512 4.807608 4.807608 4.778258 4.778258
57 58 59 60 61 62 63 64 65 66 67 68 69 70
4.474961 4.474961 4.474961 4.597350 4.597350 4.864747 4.864747 4.760193 4.586289 4.586289 4.556939 4.556939 4.253642 4.253642
71 72 73 74 75
4.253642 4.376031 4.376031 4.643428 4.643428
$se.fit
[1] 0.15659117 0.15659117 0.15232183 0.15232183 0.15150895 0.16118783 0.16118783 0.17023857 0.19749069 0.10984658 0.10984658
[12] 0.10627035 0.10627035 0.09980025 0.11627494 0.11627494 0.12476260 0.12476260 0.11171336 0.11171336 0.10903381 0.10903381
[23] 0.10000530 0.11654591 0.11654591 0.12303394 0.16348531 0.16348531 0.11170074 0.11170074 0.10884180 0.10884180 0.09841393
[34] 0.09841393 0.11378429 0.11378429 0.11938023 0.11938023 0.10769839 0.10399060 0.10399060 0.09282305 0.09282305 0.10700868
[45] 0.10700868 0.10700868 0.11275861 0.11275861 0.15407212 0.15407212 0.11205113 0.11205113 0.10690481 0.10690481 0.09678168
[56] 0.09678168 0.10691355 0.10691355 0.10691355 0.11313459 0.11313459 0.15152148 0.15152148 0.14191809 0.13638804 0.13638804
[67] 0.12964847 0.12964847 0.13426621 0.13426621 0.13426621 0.14006503 0.14006503 0.16957713 0.16957713
$df
[1] 69
$residual.scale
[1] 0.44761
diversityData$pred_shanH = predict(shanMod, se.fit = TRUE)$fit
diversityData$se_shanH = predict(shanMod, se.fit = TRUE)$se.fit
plotShannonH_pred <- diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = pred_shanH, shape = as.factor(Station), fill = as.factor(Station))) +
geom_segment(aes(y = pred_shanH - 2 * se_shanH, yend = pred_shanH + 2 * se_shanH, xend = Size_Class, color = as.factor(Station)), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"), alpha = 0.5) +
plotSpecs2 + ylab("Predicted Diversity (Shannon H)") #+ scale_y_log10()
plotShannonH_pred
predict(pielouMod, se.fit = TRUE)
$fit
1 2 3 4 5 6 7 8 9 10
0.019630086 0.019630086 0.021886750 0.021886750 0.021444103 0.024760944 0.024760944 0.022468058 0.019016295 0.010820847
11 12 13 14 15 16 17 18 19 20
0.010820847 0.013077512 0.013077512 0.012634865 0.015951706 0.015951706 0.013658820 0.013658820 0.006908498 0.006908498
21 22 23 24 25 26 27 28 29 30
0.009165163 0.009165163 0.008722516 0.012039357 0.012039357 0.009746471 0.006294707 0.006294707 0.005743555 0.005743555
31 32 33 34 35 36 37 38 39 40
0.008000220 0.008000220 0.007557573 0.007557573 0.010874414 0.010874414 0.008581527 0.008581527 0.006479321 0.008735985
41 42 43 44 45 46 47 48 49 50
0.008735985 0.008293338 0.008293338 0.011610179 0.011610179 0.011610179 0.009317293 0.009317293 0.005865530 0.005865530
51 52 53 54 55 56 57 58 59 60
0.009217729 0.009217729 0.011474393 0.011474393 0.011031746 0.011031746 0.014348587 0.014348587 0.014348587 0.012055701
61 62 63 64 65 66 67 68 69 70
0.012055701 0.008603938 0.008603938 0.013055314 0.015311978 0.015311978 0.014869331 0.014869331 0.018186172 0.018186172
71 72 73 74 75
0.018186172 0.015893286 0.015893286 0.012441523 0.012441523
$se.fit
[1] 0.005244848 0.005244848 0.005101852 0.005101852 0.005074625 0.005398808 0.005398808 0.005701953 0.006614733 0.003679190
[11] 0.003679190 0.003559408 0.003559408 0.003342699 0.003894501 0.003894501 0.004178785 0.004178785 0.003741716 0.003741716
[21] 0.003651967 0.003651967 0.003349567 0.003903576 0.003903576 0.004120886 0.005475760 0.005475760 0.003741293 0.003741293
[31] 0.003645536 0.003645536 0.003296266 0.003296266 0.003811079 0.003811079 0.003998509 0.003998509 0.003607239 0.003483050
[41] 0.003483050 0.003109006 0.003109006 0.003584138 0.003584138 0.003584138 0.003776725 0.003776725 0.005160476 0.005160476
[51] 0.003753029 0.003753029 0.003580658 0.003580658 0.003241595 0.003241595 0.003580951 0.003580951 0.003580951 0.003789318
[61] 0.003789318 0.005075045 0.005075045 0.004753390 0.004568167 0.004568167 0.004342432 0.004342432 0.004497098 0.004497098
[71] 0.004497098 0.004691323 0.004691323 0.005679799 0.005679799
$df
[1] 69
$residual.scale
[1] 0.0149922
diversityData$pred_pielouJ = predict(pielouMod, se.fit = TRUE)$fit
diversityData$se_pielouJ = predict(pielouMod, se.fit = TRUE)$se.fit
plot_pielouJ_pred <- diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = pred_pielouJ, shape = as.factor(Station), fill = as.factor(Station))) +
geom_segment(aes(y = pred_pielouJ - 2 * se_pielouJ, yend = pred_pielouJ + 2 * se_pielouJ, xend = Size_Class, color = as.factor(Station), alpha = 0.5), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both")) +
plotSpecs2 + ylab("Predicted Evenness (Pielou J)") + scale_y_log10()
plot_pielouJ_pred
plotPredictions <- plot_grid(plotObs_pred, plotChao1_pred, plotShannonH_pred, plot_pielouJ_pred, nrow = 1, labels = LETTERS)
Warning: NaNs producedWarning: Transformation introduced infinite values in continuous y-axisWarning: Removed 11 rows containing missing values (geom_segment).
plotPredictions
ggsave(here::here("Figures", "ConventionalAlphaPredictions.png"), plotPredictions, width = 11, height = 4)
plot_grid(plotObs, plotChao1, plotShan, plotPielou,
plotObs_pred, plotChao1_pred, plotShannonH_pred, plot_pielouJ_pred,
nrow = 2, labels = LETTERS)
Warning: NaNs producedWarning: Transformation introduced infinite values in continuous y-axisWarning: Removed 11 rows containing missing values (geom_segment).
alphaSummary <- tibble(
metric = c("Observed ASVs", "Richness (Chao1)", "Diversity (Shannon H)", "Evenness (Pielou J)"),
model = list(obsMod, chao1Mod, shanMod, pielouMod)
)
alphaSummary <- alphaSummary %>%
mutate(df = map(model, ~broom::tidy(summary(.))))
alphaSummary <- alphaSummary %>%
select(-model) %>%
unnest(df)
alphaSummary <- alphaSummary %>%
rename(Metric = metric, Term = term, Estimate = estimate, `Standard Error` = std.error, `T Value` = statistic, p = p.value) %>%
mutate(Term = str_replace(Term, "^I?\\((.*)\\)", "\\1"),
Term = str_replace(Term, "\\^2", "\\^2\\^"),
Term = str_replace(Term, "depth", "Depth"),
Term = str_replace(Term, "lat", "Latitude"),
Term = str_replace(Term, "_", " ")# BOOKMARK!!
) %>%
mutate(Estimate = format(Estimate, digits = 2, scientific = TRUE) %>%
reformat_sci()
) %>%
mutate(`Standard Error` = format(`Standard Error`, digits = 2, scientific = TRUE) %>%
reformat_sci()
) %>%
mutate(`T Value` = format(`T Value`, digits = 2, scientific = FALSE)) %>%
mutate(p = if_else(p < 0.001, "< 0.001", format(round(p, digits = 3)))) %>%
rename(`Standard\nError` = `Standard Error`) %>%
identity()
alphaSummary %>% flextable() %>% merge_v(j = 1) %>% theme_vanilla() %>%
bold(i = ~ p< 0.05, j = "p") %>%
colformat_md() %>%
set_table_properties(layout = "autofit") %>%
align(j = -c(1:2), align = "right")
Metric | Term | Estimate | Standard | T Value | p |
Observed ASVs | Intercept | 2.7 x 105 | 1.3 x 105 | 2.00 | 0.049 |
log(Size Class) | 3.3 x 101 | 8.1 x 100 | 4.06 | < 0.001 | |
log(Size Class)2 | -6.3 x 100 | 1.5 x 100 | -4.19 | < 0.001 | |
Latitude | -1.4 x 104 | 7.0 x 103 | -2.00 | 0.049 | |
Latitude2 | 1.8 x 102 | 9.1 x 101 | 2.00 | 0.049 | |
Depth | 4.4 x 100 | 3.2 x 100 | 1.38 | 0.173 | |
Richness (Chao1) | Intercept | 8.4 x 105 | 4.2 x 105 | 2.00 | 0.049 |
log(Size Class) | 7.8 x 101 | 2.5 x 101 | 3.07 | 0.003 | |
log(Size Class)2 | -1.6 x 101 | 4.7 x 100 | -3.38 | 0.001 | |
Latitude | -4.4 x 104 | 2.2 x 104 | -2.00 | 0.049 | |
Latitude2 | 5.7 x 102 | 2.9 x 102 | 2.00 | 0.049 | |
Depth | 1.8 x 101 | 1.0 x 101 | 1.80 | 0.077 | |
Diversity (Shannon H) | Intercept | 1.7 x 103 | 7.9 x 102 | 2.12 | 0.037 |
log(Size Class) | 2.0 x 10-1 | 4.8 x 10-2 | 4.25 | < 0.001 | |
log(Size Class)2 | -3.7 x 10-2 | 8.9 x 10-3 | -4.15 | < 0.001 | |
Latitude | -8.7 x 101 | 4.1 x 101 | -2.11 | 0.038 | |
Latitude2 | 1.1 x 100 | 5.3 x 10-1 | 2.11 | 0.038 | |
Depth | 1.9 x 10-2 | 1.9 x 10-2 | 0.99 | 0.325 | |
Evenness (Pielou J) | Intercept | -2.2 x 101 | 2.6 x 101 | -0.84 | 0.405 |
log(Size Class) | -4.0 x 10-3 | 1.6 x 10-3 | -2.48 | 0.016 | |
log(Size Class)2 | 6.8 x 10-4 | 3.0 x 10-4 | 2.28 | 0.026 | |
Latitude | 1.1 x 100 | 1.4 x 100 | 0.84 | 0.405 | |
Latitude2 | -1.5 x 10-2 | 1.8 x 10-2 | -0.84 | 0.406 | |
Depth | -3.0 x 10-4 | 6.3 x 10-4 | -0.48 | 0.636 |
richSummary %>% rename_(.dots = setNames(names(.), paste0('break_', names(.))))
Warning: `rename_()` was deprecated in dplyr 0.7.0.
Please use `rename()` instead.
diversityDataWB <- full_join(diversityData,
richSummary %>% rename_(.dots = setNames(names(.), paste0('break_', names(.)))),
by = c("ID" = "break_sample_names"), suffix = c("", "_break")) %>%
mutate(pielouJ2 = shannonH/break_estimate) %>%
identity()
diversityDataWB
pielouMod2 <- lm(pielouJ2 ~ log(Size_Class) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityDataWB)
summary(pielouMod2)
Call:
lm(formula = pielouJ2 ~ log(Size_Class) + I(log(Size_Class)^2) +
lat + I(lat^2) + depth, data = diversityDataWB)
Residuals:
Min 1Q Median 3Q Max
-0.013930 -0.005148 -0.002493 0.000898 0.105967
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.957e+01 2.755e+01 -0.710 0.4799
log(Size_Class) -3.287e-03 1.669e-03 -1.969 0.0530 .
I(log(Size_Class)^2) 5.741e-04 3.106e-04 1.848 0.0688 .
lat 1.018e+00 1.436e+00 0.709 0.4808
I(lat^2) -1.322e-02 1.869e-02 -0.707 0.4817
depth -2.380e-04 6.612e-04 -0.360 0.7200
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.01569 on 69 degrees of freedom
Multiple R-squared: 0.06696, Adjusted R-squared: -0.000652
F-statistic: 0.9904 on 5 and 69 DF, p-value: 0.4301
Ok. So the narrative makes sense. Alpha diveristy is driven by variability in richness rather than evenness. Why would we see an effect in chao1 but not breakaway? Because chao1 is more driven by abundant stuff that makes the rarification threshold. My first hunch is that chao1 responds to evenness, but actually that shouldn’t have any effect since there is no evenness variability? Or maybe just that breakaway is more variable (because it detects fine level differences in rare species) and that doesn’t map as nicely with overall patterns.
plotBreak <- diversityDataWB %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = break_estimate, shape = as.factor(Station), fill = as.factor(Station))) +
plotSpecs +
#scale_y_log10()+
ylab("Richness (Breakaway)")
plotBreak
plotPielou2 <- diversityDataWB %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = pielouJ2, shape = as.factor(Station), fill = as.factor(Station))) +
plotSpecs +
#scale_y_log10()+
ylab("Evenness (PielouJ)")
plotPielou2
predict(pielouMod2, se.fit = TRUE)
$fit
1 2 3 4 5 6 7 8 9
0.0116843027 0.0116843027 0.0135766875 0.0135766875 0.0133698701 0.0160079436 0.0160079436 0.0139748525 0.0099214190
10 11 12 13 14 15 16 17 18
0.0043268153 0.0043268153 0.0062192001 0.0062192001 0.0060123827 0.0086504562 0.0086504562 0.0066173651 0.0066173651
19 20 21 22 23 24 25 26 27
0.0011039449 0.0011039449 0.0029963297 0.0029963297 0.0027895123 0.0054275858 0.0054275858 0.0033944947 -0.0006589388
28 29 30 31 32 33 34 35 36
-0.0006589388 0.0002124081 0.0002124081 0.0021047929 0.0021047929 0.0018979755 0.0018979755 0.0045360490 0.0045360490
37 38 39 40 41 42 43 44 45
0.0025029579 0.0025029579 0.0009065649 0.0027989497 0.0027989497 0.0025921323 0.0025921323 0.0052302058 0.0052302058
46 47 48 49 50 51 52 53 54
0.0052302058 0.0031971147 0.0031971147 -0.0008563188 -0.0008563188 0.0033197651 0.0033197651 0.0052121500 0.0052121500
55 56 57 58 59 60 61 62 63
0.0050053325 0.0050053325 0.0076434061 0.0076434061 0.0076434061 0.0056103150 0.0056103150 0.0015568815 0.0015568815
64 65 66 67 68 69 70 71 72
0.0066525974 0.0085449823 0.0085449823 0.0083381648 0.0083381648 0.0109762384 0.0109762384 0.0109762384 0.0089431473
73 74 75
0.0089431473 0.0048897138 0.0048897138
$se.fit
[1] 0.005490645 0.005490645 0.005340947 0.005340947 0.005312444 0.005651821 0.005651821 0.005969172 0.006924728 0.003851613
[11] 0.003851613 0.003726218 0.003726218 0.003499353 0.004077014 0.004077014 0.004374622 0.004374622 0.003917069 0.003917069
[21] 0.003823115 0.003823115 0.003506543 0.004086515 0.004086515 0.004314009 0.005732378 0.005732378 0.003916626 0.003916626
[31] 0.003816382 0.003816382 0.003450743 0.003450743 0.003989683 0.003989683 0.004185897 0.004185897 0.003776290 0.003646281
[41] 0.003646281 0.003254708 0.003254708 0.003752106 0.003752106 0.003752106 0.003953719 0.003953719 0.005402318 0.005402318
[51] 0.003928913 0.003928913 0.003748464 0.003748464 0.003393511 0.003393511 0.003748771 0.003748771 0.003748771 0.003966902
[61] 0.003966902 0.005312884 0.005312884 0.004976155 0.004782251 0.004782251 0.004545938 0.004545938 0.004707852 0.004707852
[71] 0.004707852 0.004911180 0.004911180 0.005945979 0.005945979
$df
[1] 69
$residual.scale
[1] 0.0156948
diversityDataWB$pred_pielouJ2 = predict(pielouMod2, se.fit = TRUE)$fit
diversityDataWB$se_pielouJ2 = predict(pielouMod2, se.fit = TRUE)$se.fit
plot_pielouJ2_pred <- diversityDataWB %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = pred_pielouJ2, shape = as.factor(Station), fill = as.factor(Station))) +
geom_segment(aes(y = pred_pielouJ2 - 2 * se_pielouJ2, yend = pred_pielouJ2 + 2 * se_pielouJ2, xend = Size_Class, color = as.factor(Station), alpha = 0.5), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both")) +
plotSpecs2 + ylab("Predicted Evenness (Pielou J2)") #+ scale_y_log10()
plot_pielouJ2_pred
plotBreakaway <- diversityDataWB %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = break_estimate, shape = as.factor(Station), fill = as.factor(Station))) +
plotSpecs +
geom_errorbar(aes(ymin = break_lower, ymax = break_upper), width = -.1) +
scale_y_log10() +
ylab("Richness (Breakaway)")
plotBreakaway
#predict(breakLm, se.fit = TRUE)
# doesn't work because built with a different data frame
Why are these not smooth curves?!! What if I redo the model, this time with the same data frame
breakLm2 <- lm(break_estimate ~ log(Size_Class) + I(log(Size_Class) ^2) + lat + I(lat^2) + depth ,data = diversityDataWB)
breakLm2 %>% summary()
Call:
lm(formula = break_estimate ~ log(Size_Class) + I(log(Size_Class)^2) +
lat + I(lat^2) + depth, data = diversityDataWB)
Residuals:
Min 1Q Median 3Q Max
-2974.5 -1191.2 -151.6 599.9 6768.1
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7124615.61 3339862.88 2.133 0.0365 *
log(Size_Class) 244.45 202.35 1.208 0.2312
I(log(Size_Class)^2) -75.16 37.65 -1.996 0.0498 *
lat -370568.38 174035.93 -2.129 0.0368 *
I(lat^2) 4817.28 2265.61 2.126 0.0371 *
depth 151.10 80.15 1.885 0.0636 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1902 on 69 degrees of freedom
Multiple R-squared: 0.1414, Adjusted R-squared: 0.0792
F-statistic: 2.273 on 5 and 69 DF, p-value: 0.0567
Note the non statistical significance overall
#predict(breakLm2, se.fit = TRUE)
diversityDataWB$pred_break = predict(breakLm2, se.fit = TRUE)$fit
diversityDataWB$se_break = predict(breakLm2, se.fit = TRUE)$se.fit
plot_break_pred <- diversityDataWB %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
# filter(Station == 4.3) %>%
ggplot(aes(x = Size_Class, y = pred_break, shape = as.factor(Station), fill = as.factor(Station))) +
geom_segment(aes(y = pred_break - 2 * se_break, yend = pred_break + 2 * se_break, xend = Size_Class, color = as.factor(Station), alpha = 0.5), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both")) +
plotSpecs2 + ylab("Predicted Richness (Breakaway -- LM)") #+ scale_y_log10()
plot_break_pred
plotAlphaWB <- plot_grid(plotBreakaway, plotShan, plotPielou2, nrow = 1, labels = LETTERS)
plotAlphaWB
ggsave(here::here("Figures", "BreakawayAlpha.png"), plotAlpha, width = 11, height = 4)
Summary table I want both breakaway metrics here
bettaTable <- myBet$table %>%
as.data.frame() %>%
rename(estimate = Estimates,
`std.error` = `Standard Errors`,
`p.value`=`p-values`
) %>%
mutate(`statistic` = NA) %>%
rownames_to_column(var = "term") %>%
select(term, estimate, std.error, statistic, p.value) %>%
as_tibble()
bettaTable
alphaSummary2 <- tibble(
metric = c("Richness (Breakaway -- LM)", "Diversity (Shannon H)", "Evenness (Pielou J)"),
model = list(breakLm, shanMod, pielouMod2)
)
alphaSummary2 <- alphaSummary2 %>%
mutate(df = map(model, ~broom::tidy(summary(.))))
## Add in willis variables
breakawaySummary <- tibble(
metric = "Richness (Breakaway -- Betta)",
model = NULL,
df = list(bettaTable)
)
alphaSummary2 = bind_rows(breakawaySummary, alphaSummary2)
alphaSummary2 <- alphaSummary2 %>%
select(-model) %>%
unnest(df)
alphaSummary2 <- alphaSummary2 %>%
rename(Metric = metric, Term = term, Estimate = estimate, `Standard Error` = std.error, `T Value` = statistic, p = p.value) %>%
mutate(Term = str_replace(Term, "^I?\\((.*)\\)", "\\1"),
Term = str_replace(Term, "\\^2", "\\^2\\^"),
Term = str_replace(Term, "depth", "Depth"),
Term = str_replace(Term, "lat", "Latitude"),
Term = str_replace(Term, "_", " ")# BOOKMARK!!
) %>%
mutate(Estimate = format(Estimate, digits = 2, scientific = TRUE) %>%
reformat_sci()
) %>%
mutate(`Standard Error` = format(`Standard Error`, digits = 2, scientific = TRUE) %>%
reformat_sci()
) %>%
mutate(`T Value` = format(`T Value`, digits = 2, scientific = FALSE)) %>%
mutate(p = if_else(p < 0.001, "< 0.001", format(round(p, digits = 3)))) %>%
rename(`Standard\nError` = `Standard Error`) %>%
identity()
alphaSummary2
alphaTable2 <- alphaSummary2 %>% flextable() %>% merge_v(j = 1) %>% theme_vanilla() %>% bold(i = ~ p< 0.05, j = "p") %>% colformat_md() %>% set_table_properties(layout = "autofit") %>%
align(j = -c(1:2), align = "right")
alphaTable2
Metric | Term | Estimate | Standard | T Value | p |
Richness (Breakaway – Betta) | Intercept | 7.1 x 106 | 2.4 x 102 | NA | < 0.001 |
log(Size Class) | 1.2 x 102 | 6.1 x 101 | NA | 0.058 | |
log(Size Class)2 | -5.0 x 101 | 1.2 x 101 | NA | < 0.001 | |
Latitude | -3.7 x 105 | 6.1 x 100 | NA | < 0.001 | |
Latitude2 | 4.8 x 103 | 1.6 x 10-1 | NA | < 0.001 | |
Depth | 1.5 x 102 | 1.0 x 101 | NA | < 0.001 | |
Richness (Breakaway – LM) | Intercept | 7.1 x 106 | 3.3 x 106 | 2.13 | 0.036 |
log(Size Class) | 2.4 x 102 | 2.0 x 102 | 1.21 | 0.231 | |
log(Size Class)2 | -7.5 x 101 | 3.8 x 101 | -2.00 | 0.050 | |
Latitude | -3.7 x 105 | 1.7 x 105 | -2.13 | 0.037 | |
Latitude2 | 4.8 x 103 | 2.3 x 103 | 2.13 | 0.037 | |
Depth | 1.5 x 102 | 8.0 x 101 | 1.89 | 0.064 | |
Diversity (Shannon H) | Intercept | 1.7 x 103 | 7.9 x 102 | 2.12 | 0.037 |
log(Size Class) | 2.0 x 10-1 | 4.8 x 10-2 | 4.25 | < 0.001 | |
log(Size Class)2 | -3.7 x 10-2 | 8.9 x 10-3 | -4.15 | < 0.001 | |
Latitude | -8.7 x 101 | 4.1 x 101 | -2.11 | 0.038 | |
Latitude2 | 1.1 x 100 | 5.3 x 10-1 | 2.11 | 0.038 | |
Depth | 1.9 x 10-2 | 1.9 x 10-2 | 0.99 | 0.325 | |
Evenness (Pielou J) | Intercept | -2.0 x 101 | 2.8 x 101 | -0.71 | 0.480 |
log(Size Class) | -3.3 x 10-3 | 1.7 x 10-3 | -1.97 | 0.053 | |
log(Size Class)2 | 5.7 x 10-4 | 3.1 x 10-4 | 1.85 | 0.069 | |
Latitude | 1.0 x 100 | 1.4 x 100 | 0.71 | 0.481 | |
Latitude2 | -1.3 x 10-2 | 1.9 x 10-2 | -0.71 | 0.482 | |
Depth | -2.4 x 10-4 | 6.6 x 10-4 | -0.36 | 0.720 |
alphaTable2 %>% save_as_docx(path = here::here("Tables", "alphaTable2.docx"))
myBet$table
plot_grid(plot_break_pred,plotShannonH_pred,plot_pielouJ2_pred, nrow = 1, labels = LETTERS)